# TODO: Add comment
# 
# Author: guochun
###############################################################################

source("./sourceAll.R")
data=getData("bci5.txt")
availible=data$status=="A" & !is.na(data$gx) & !is.na(data$gy) &
		data$gx <50 & data$gy < 50 
subdata=data[availible,c(1,3:5)]
subdata$dbh=subdata$dbh/1000

#change the single tree point to 32 point around the tree center
alpha=seq(0,2*pi,pi/16)
alpha=alpha[-length(alpha)]

newx=as.vector(apply(subdata,1,function(x) x[2]+x[4]*cos(alpha)))
newy=as.vector(apply(subdata,1,function(x) x[3]+x[4]*sin(alpha)))
newtag=rep(subdata$tag,each=length(alpha))

#选择一个站点坐标
totalsitexy=expand.grid(seq(0,50,length=10),seq(0,50,length=10))
totalsampled=rep(FALSE,length(newx))

for(k in 1:dim(totalsitexy)[1]){
	temp=lightSample(totalsitexy[k,1],totalsitexy[k,2],newx,newy,subdata,30)
	totalsampled=totalsampled + temp$sampled
	print(k)
}

sum(totalsampled)/length(totalsampled)











lightSample=function(x.site,y.site,newx,newy,subdata,scanr){

#坐标原点转换为站中心
newx=newx-x.site
newy=newy-y.site
subdata$gx=subdata$gx-x.site
subdata$gy=subdata$gy-y.site
#将xy坐标转换为以基站为原点的极坐标
r=sqrt(newx^2+newy^2)
beta=asin(newy/r)
newdata=data.frame(tag=newtag,r=r,beta=beta,sampled=rep(FALSE,time=length(beta)))
r0=sqrt(subdata$gx^2+subdata$gy^2)
beta0=asin(subdata$gy/r0)
newdata0=data.frame(tag=subdata$tag,r0=r0,beta0=beta0,sampled=rep(FALSE,time=length(r)))
#由于有最大扫描半径，所以要先设置需要扫描的个体
nmax=sum(r0<scanr)

newdata0=newdata0[order(newdata0$r0),]
#保留排序后与原来顺序的关系
#reorder=match(r0,newdata0$r0)
#全局遮挡条件
cond=data.frame(startAngle=NA,endAngle=NA,maxr=NA)

for(i in 1:nmax){
	#从最近开始选择一棵树
	tree=newdata0[i,]
	nwhich=which(newdata$tag==tree$tag)
	treepoints=newdata[nwhich,]
	#计算可视范围内的点
	if(!is.na(cond[1,1])){
		for(j in 1:dim(cond)[1]){
			pre.cond=cond[j,]
			unsee=treepoints$r > pre.cond$maxr & 
					treepoints$beta > pre.cond$startAngle & 
					treepoints$beta < pre.cond$endAngle
			if(any(unsee)){
				treepoints=treepoints[-which(unsee),]
				nwhich=nwhich[-which(unsee)]
			}
				
		}
			
	}
	#当还有未遮挡点时，计算新的遮挡区域和照射到的点
	if(dim(treepoints)[1]>0){
		#计算遮挡条件
		cond.angle=range(treepoints$beta)
		cond.r=max(treepoints$r[treepoints$beta==cond.angle[1] | treepoints$beta==cond.angle[2]])
		#储存新的遮挡条件
		if(is.na(cond[1,1])){
			cond[1,]=c(cond.angle,cond.r)
		}
		else{
			cond=rbind(cond,c(cond.angle,cond.r))
		}
		#将照射到的点设成sampled=T
		sampled=nwhich[treepoints$r<=cond.r & treepoints$beta >= cond.angle[1] & treepoints$beta <= cond.angle[2]]
		newdata$sampled[sampled]=T
	}
	
}
#newdata=newdata[reorder,]
return(newdata)
}
